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Abstract 

We study the phase diagram of the four dimensional Ising model with first 
and second neighbour couplings, specially in the antiferromagnetic region, by 
using Mean Field and Monte Carlo methods. From the later, all the transition 
lines seem to be first order except that between ferromagnetic and disordered 
phases in a region including the first-neighbour Ising transition point. 
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Antiferromagnetism has been considered in a great variety of models in order 
to find properties not present in the purely ferromagnetic (FM) system. Mostly 
they have been models in two and three dimensions. For instance, in several works 
on High Temperature Superconductivity the transition from paramagnetic (PM) to 
non-pure FM ordered phases has been studied [I], 0, [5[ On the other hand, in 
relation to the finite temperature phase transition in pure gauge 577(3), the d = 3 
three state Potts model with negative second neighbour couplings was unsuccessfully 
considered to find a new critical behaviour Q . 

In four dimensions, in diluted systems recently new critical exponents have been 
obtained 0. From a quantum field theory point of view, no argument appears to 
prevent the existence of a non-trivial ultraviolet limit in an antiferromagnetic (AF) 
lattice 4 theory |Q. In fact, Gallavotti and Rivasseau more than ten years ago 
considered the possibility of an AF action which could change the ultraviolet limit 
of the pure 4 model ||. 

Motivated in part by that, we decided to study the phase diagram and possible 
critical behaviour of a d — 4 Ising model with an AF phase non trivially equivalent 
to the standard FM one. As we will see, we have found a rich phase diagram with 
several phase transitions but none of them seems to show a new critical behaviour. 

The naive way to introduce antiferromagnetism in the Ising model is to consider 
a negative coupling. In this case the state with minimal energy for large (3 is a 
staggered vacuum. On a hypercubic lattice, we denote the coordinates of site n as 
(n x ,n y ,n z ,n t ). If we make the transformation 

a(n x , n y , n z , n t ) -> (-l) n * +n « +n * +nt a(n x , n y , n z , n t ), (1) 

the system with negative (3 is mapped onto the positive (3 one, both regions being 
exactly equivalent. 

Therefore, to consider true antiferromagnetism we must take into account either 
different geometries or more couplings, in order to make the transformation (||) not 
an exact mapping. The option we have chosen here is to add a coupling between 
points at a distance of v2 lattice units. 

So, we will work with an Ising model with Hamiltonian 

H = /3i^a(n)a(n + /2) + / 9 2 J2 o{n)a{n + p + v), (2) 

on a four dimensional hypercubic lattice, with side L and periodic boundary condi- 
tions. Here jl denotes the unitary vector in the \x direction. 

The transformation ([]]) maps the semiplane f3\ > onto the f3± < one, and 
therefore only the region with j3i > will be considered. On the line /3i — the 
system decouples into two independent sublattices. 

The presence of two couplings with opposite signs makes frustration to appear, 
and very different vacua are possible. For small values of /3i and /3 2 the system is 
disordered (PM phase). On the other hand, we have computed the configurations 
which minimize the energy for several asymptotic values of the parameters. We have 
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only considered configurations with periodicity two. More complex structures have 
not been observed in our simulations. 
We have found the following regions: 

1. Small absolute values of /3 2 - PM phase. 

2. Large positive /3i, /3 2 , or large f3\ > 0, (3 2 < with f3\ > 6\p 2 \. The vacuum is 
the FM one, a(n) = a (a stands for a fixed spin). 

3. Large fix > 0,{3 2 < with fi 1 < 6\f3 2 \ and fi 1 > 2\fi 2 \. 111 this region the 
vacuum is an FM configuration on a three-dimensional cube and AF on the 
other direction ii. We have 4 identical possibilities to break the symmetry: one 
for every possible choice of ii (more precisely 8 if we consider also the global 
cr(n) — > —cr{n) symmetry). We call this vacuum Hyperplane Antiferromag- 
netic (HPAF). The configuration is of type a(n) = (— 1)™ M (T , where \x can be 
any direction. 

4. Large f3\ > 0, (3 2 < with i3\ < 2\(3 2 \. In this region the vacuum is an FM con- 
figuration on a two-dimensional plane and AF on the other two directions. We 
have now six equivalent vacua. We call this vacuum Plane Antiferromagnetic 
(PAF), and the configuration can be written as cr(n) = (— l) rifl+ri,/ o"o, where 
fi, v (fi < v) can be any direction. 

We remark that, in order to avoid undesirable (frustrating) boundary effects for 
ordered phases, we work with even lattice side L as periodic boundary conditions 
are imposed. 

Now we must define an order parameter for every phase. 

For the FM phase the order parameter is the standard magnetization 



where V is the lattice volume. 

In the HPAF phase we define the parameter M 2 ^ = (l/V) En(-1)™ M °"(^)- M 2 ^ 
will be different from zero only in the HPAF phase, where the system becomes 
antiferromagnetic on the ji direction. We have 4 order parameters (one for every 
possible value of li) and only one of them will be different from zero in the HPAF 
phase. To consider only one, we define 



that, although it is not a true order parameter, is more appropriate for measuring 
on a finite lattice. 



For the PAF region we first define M 3j ^ u = (l/V) En(-l) n " +n ''o'(n)' The same 
previous discussion applies, and we will work with 



(3) 




(4) 




(5) 
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To understand the behaviour of the system, we have carried out a Mean Field 
analysis and a Monte Carlo simulation. 

From the Mean Field analysis we can plot a general phase diagram. We use 
the standard technique, also used for gauge theories, fllP)] . We define three different 
order parameters Vy, V 2 , V 3 , labelling the FM, HPAF and PAF phases respectively. 

Our Mean Field ansatz is a combination of the three possible order parameters 

V(n) = Vy + {-l) nt V 2 + {-l) n * +n yV 3 , (6) 

with their corresponding auxiliary fields A(n) = Ai + (— l) nt v42 + (— l) nx+ny A 3 , where 
we have selected a fixed breaking direction for the HPAF and the PAF phases. 
After a standard computation, the free energy per site becomes 



-[(4A + 12/3 2 )V? + 2ft V 2 2 - A(3 2 V£ 
+A 1 V 1 + A 2 V 2 + A 3 V 3 + T], 



(7) 



with T defined as 



(8) 



T = (1/4) [log cosh(A! + A 2 + A 3 ) + log cosh(Ai + A 2 - A 3 ) 
+ logcosh(A 1 - A 2 + A 3 ) + logcosh^! - A 2 - A 3 )}. 

By solving the saddle point equations for F 
dF dF 

wr ' air - i=1 ' 2 ' 3 ' (9) 

we will find regions in the parameter space where the minimum of F corresponds to 
different values of V^. 

This system of 6 equations is easily reduced to 3 equations resembling the clas- 
sical M = a tanh(M) but they are coupled. By inspection one finds some solutions: 

1. Vy = V 2 = V3 = is always a solution. 

2. Vy ^ 0, V 2 = V 3 = is solution if 8ft + 24ft > 1. 

3. V 2 ^ 0, Vy = V 3 = if 4ft > 1. 

4. V 3 =£0,Vy = V 2 = if -8ft > 1. 



This divides the parameter space in several regions, some of them being mixed. 
In every region one can see that those are the unique solutions, and the deeper 
minimum gives the true solution for each pair (ft, ft). The minimum with V. — Vi 
corresponds to the PM phase. If the minimum is at Vy ^ 0, V 2 — V 3 — 0, we are 
in the FM phase. If only V 2 7^ 0, this corresponds to the HPAF phase and, if only 
V 3 7^ 0, the phase is PAF. The result for the phase diagram is shown in Fig. 1. 

The transitions separating the different phases are straight lines: 8ft + 24ft = 1 
between PM and FM, 4ft = 1 for PM-HPAF, -8ft = 1 for PM-PAF, ft + 6ft = 
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for FM-HPAF and f3 1 + 2f3 2 = for PAF-HPAF. The lines between any ordered 
phase and the PM one are second order with the classical Mean Field exponents 
(the shape of the equations defining these lines is identical for all of them). The 
lines between two ordered phases are first order (when crossing the lines an order 
parameter jumps abruptly from to a positive value). 

After this Mean Field approach, we have run a Heat Bath Monte Carlo compu- 
tation. 

We measure the energies, defined as 

1 dlogl 1 dlogZ 

{El) ~ I7^T' {E2) ~ l27^T' (10) 

and the expectation value of the order parameters, (Mi), (M 2 ), (M 3 ). 

With these quantities we have studied the global phase diagram. Using the 



spectral density method [JTTJ] and hysteresis cycles on lattices of size L = 8, we have 
found the transition lines shown in Fig. 1. 

The line FM-PM which includes the standard Ising model point (/3 2 = 0) is 
second order with mean field exponents (from fli — to some f3i larger than the 
cited standard critical point). 

The transition lines FM-HPAF and HPAF-PAF (see Fig. 1) for large values of j3\ 
behave as strong first order transitions: they present metastables states and evolve 
very slowly with our local Monte Carlo simulations. Concerning this, we address 
the question whether the ordered phases are directly connected or there exists a PM 
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region between them. We could expect a PM region separating the ordered phases 
along the asymptotes Pi + 6/3 2 = and Pi + 2/3 2 = (narrower for larger Pi values) 
because the ground state is highly degenerate for those values of Pi, P2 (you can place 
an FM sublattice beside an HPAF one without the interfase increasing the energy, 
if you choose the correct orientation. You can also do this for the HPAF-PAF case). 
However, from our MC simulation it is not possible to give a conclusive answer since 
the width of the hypothetical PM region decreases when increasing Pi, and for a 
fixed lattice size there is a practical limit in the precision of the measures of critical 
values. We have found that the lines PM-FM and PM-HPAF approach very fast 
and they cannot be resolved for large Pi values. A similar situation is found when 
the PM-HPAF line comes near the PM-PAF one. In this case the approaching is 
even faster. 

The transition PM-HPAF presents a clear metastability with a large latent heat, 
indicating a first order transition. 

The transition PM-PAF has a very involved behaviour. For large values of Pi 
the transition shows a large latent heat. When we move towards smaller Pi this 
latent heat decreases, and for values Pi < 0.20 it disappears when measured from 
hysteresis cycles on lattices L = 8, 12, 16. This led us to think about the possibility 
of having a second order transition and study it in a more detailed way. 

On this line, we have studied three points: Pi = 0.1, Pi = 0.05 and Pi = 0.0. 

At Pi = 0.0 one has two decoupled lattices (one being the constituted by all the 
first neighbours of the other one). Then we have simulated separately these lattices 
(we have run on FA lattices). In the following, when we talk about the size of the 
lattice L at the point Pi = 0.0, we will be describing an FA lattice with a number 
of sites L A /2. 

We have simulated lattices from L = 6 to L = 24. The number of Monte Carlo 
iterations has depended on the lattice size, going approximately from two hundred 
thousand (for L = 6) to more than one million (for L = 20, 24), discarding a quantity 
of the order of a twenty per cent for thermalization. The largest autocorrelation time 
found has been O(10 3 ) for large lattices. 

We have used the spectral density method to locate the transition points. No 
signal of metastability has appeared for the smallest lattices. However on the largest 
lattices one can observe some trace of a two peak structure. At Pi = 0.1 one needs 
an L = 20 lattice to begin to distinguish those signs and at Pi = 0.05 one has to 
go up to L = 24. However at Pi = 0.0 one can already see a two peak structure at 
L = 16. 

In Fig. 2 we show energy histograms at the three points for different lattice sizes. 

Besides this direct look at the histograms, we have carried out at those points a 
Finite Size Scaling analysis [|12|, [13j using lattices from L = 6 to L = 24 in order to 
clarify the order of the transition. Despite one could not strictly talk about critical 
exponents if the transition resulted to be first order, it is always possible to calculate 
the exponents with which some quantities diverge with increasing L and see if they 
agree with those expected for a first order transition {y = 1/d = 0.25, a = 1, 7 = 1). 
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L=16 L=20 L=24 




L=12 L=16 L=20 L=24 




Figure 2: Energy histograms for different lattice sizes in the proximities of the three 
transition points studied on the line PM-PAF. 



We have computed the exponents 7 and v in the way exposed below. A good 
estimation of a has not been possible. We have used (Ei), (£2), the relevant order 
parameter which is the expectation value of (M 3 ) ((M) hereafter) and some of their 
derivatives. 

We start obtaining v from the scaling law of the quantity 



which is 



, max 



(L) ~ L 1/u . (12) 
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To obtain 7 we use the susceptibility 

X = V((M 2 } — (M) 2 ), (13) 

with scaling law 

Xmax (L) ~ L?l v . (14) 

The behaviour of the maximal value of k 2 with L, from where we extract u, is 
shown in Fig. 3. 

From this figure, it is clear that not all the L values simulated are inside the 
asymptotic region where one could define a correct v exponent. Contrarily, one ob- 
serves that the exponent varies if obtained from different sets of lattice sizes. Notice 
that the behaviour for L < 20 is surprisingly linear, with a slope corresponding to 
v > 1/d. However the L = 20, 24 lattices point to a value v ~ l/d. 

A similar behaviour is found in the computation of 7/V. 

The results for the exponents obtained separately from small and large lattices 
are shown in Table 1. 

On the other hand, we tried to obtain a/u from the scaling of the elements of 
the matrix 

C»=m (15) 

Its eigenvectors should be respectively orthogonal and parallel to the direction of 
the transition line. This line is, in this region, almost parallel to the f3\ axis and 
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/3i = 0.0 


Px = 0.05 


Pi = 0.1 




1 st ord 


L=8,12,16 


L=20,24 


L=8,12,16 


L=20,24 


L=6,8,10,16 


L=20,24 


V 


0.25 


0.359(4) 


0.245(33) 


0.346(7) 


0.223(16) 


0.342(6) 


0.224(14) 


1 


1 


1.08(2) 


1.06(23) 


1.08(3) 


1.03(10) 


1.08(3) 


1.02(9) 



Table 1: Exponents obtained at the three transition points studied on the line PM- 
PAF using different sets of lattice sizes. The first column refers to the expected for 
a first order transition. 



V = 1 -^h- (16) 



the C 2,2 element results to be very close to the eigenvalue corresponding to the 
orthogonal eigenvector, which will be the only relevant element of the matrix C M (it 
would not be the case if the point were a multicritical point). In short, we studied 
the divergence of the C 2,2 matrix element. However we did not get a reasonable 
estimation of a. There are basically two coupled difficulties. Firstly the divergence 
of the specific heat is usually difficult to fit. Many times one does not need just a 
power term in L but also a constant term or even a logarithmic or exponential one. 
This problem could be solved but the added objection here is that we have very 
little data inside the asymptotic region, and the estimations we can make for a are 
too much poor. 

A satisfactory conclusion is obtained from the behaviour of the Binder cumulant 

M 

(El) 

my 

The minimum of this observable approaches the value 2/3 with increasing L in a 
second order transition. This is not the case if the transition is first order, where 
in the thermodynamical limit it takes a value depending on the position of the two 
energy peaks and easily determined from ([16]). If these peaks are delta functions 
situated at E a and E b , then 

2{E a i + E b 4 ) 

► mm — 1 ,n o, 2 • \ x 1 / 

3(E a 2 + E b 2 f 

In Fig. 4 it is shown how, despite its closeness to 2/3, V min tends with increasing 
L to a value slightly but significantly smaller than this (approximately 0.6658) which 
indicates a first order transition with a latent heat that can be estimated in the 
following way. Let write E b = E a + A. Looking at the histograms (Fig. 2) one 
observes that, approximately, one peak is situated at an energy of —0.1 and the 
latent heat A is one order of magnitude smaller. Taking that value for E a and the 
value 0.6658 for V m i n in the thermodynamic limit, one obtains from (|T7| ) a latent 
heat A ~ —0.005 which is of the same order as the expected from the histograms. 

On the other hand, we have computed the transition points /3|(oo) at each (5\ 
value by using the expression 

/5 2 c (oo) = (3 C 2 (L) - AL- l ' v . (18) 
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0.667 




1/L 4 

Figure 4: Minimum of the Binder cumulant as a function of the lattice size. The 
cross marks the value 2/3. 



We have used only the largest values of L and the corresponding exponent v. 
For each observable, we have (in principle) a different /3|(L) obtained as the po- 
sition of the maximum derivative (with respect to fa). The resulting /3f(oo) is 
very similar for any of them and we take the mean value. We obtain /3f(cx>) = 
-0.17904(3), -0.17568(8), -0.17459(15) at fa = 0.1,0.05,0.0 respectively. 

As it was met in phase transitions of some other systems with frustration effects 
and non-trivial sublattice structures, for instance [f|, [HJ, this study has revealed to 



be hard and we have needed rather large lattices to see the asymptotic behaviour, 
which was masked up to a significant size (L = 24 is a considerably large lattice in 
d = 4). On the other hand, a usual problem in this kind of systems is that they are 
difficult to equilibrate. However, we have simulated 64 independent lattices starting 
from different configurations and, after thermalization, they give fully compatible 
results. 

In conclusion, we have described the phase diagram of the four dimensional Ising 
model with first, fa, and second, fa, neighbour couplings. In the fa > semiplane, 
four regions (PM, FM, HPAF and PAF) are present. The fa < semiplane is 
obtained from the above by means of a trivial transformation. All the transition 
lines have resulted to be first order except that corresponding to the first neighbour 
Ising model which is second order with Mean Field exponents. In consequence, this 
simple modification to the Ising model does not seem, by itself, to be useful for 
finding a critical behaviour different from the classical one. 
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